%read an display global ocean mag data
%29.11.2007
%fid = fopen('C:\Manoj\ocean\mf6\1x1_ecco_mean_2000_2006.res','rb')
fid = fopen('C:\Manoj\ocean\agu07\ecco_ann_equinox.source.res','rb')
%fid = fopen('C:\Manoj\ocean\temp\ecco_annual_only_imag1.res','rb')
fread(fid,518400,'uchar'); %Header includes the gird size, files used for computation etc
[Hpr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi; %in nT
[Hpi(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Htr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hti(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hrr(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
[Hri(:,:)] = reshape(fread(fid,360*180,'float32'),[180,360])*400*pi;
fclose all

load coast;
data = flipud(Hti);
worldmap('world');
h=contourfm(data,[1,89.9,0.5],20);
geoshow(lat, long,'color','k');
caxis([-.05 .05]);
title('Im B_{\theta}');
colorbar;
%map=colormap('gray');
%colormap(flipud(map));
% %caxis([-.1,.1]);
% title('Im Br ECCO 2000-2006 Annual at 400 km');
% colorbar('horiz');
% %title('Annual harmonics of Meridional V component of ECCO 2000-20006');
% %title('Bi-annual (6 months) harmonics of Meridional V component of ECCO 2000-20006');
% title('monthly harmonics of Meridional V component of ECCO 2000-20006');
% 
% %write the results to a file
% 
load c:\manoj\ocean\ECCO_Cord;
[X1,Y1] = meshgrid(V_long,[89.5:-1:80.5 fliplr(U_lat') -80.5:-1:-89.5]);
 
fid = fopen('c:\manoj\ocean\agu07\ECCO_2000_2006_ANNUAL_B.txt','wt');
fprintf(fid,'Lat   Long    Hrr     Hri     Htr     Hti     Hpr     Hpi\n');
% fprintf(fid,'Lat   Long    Hrr     Htr     Hpr \n');
 for i = 1:360,
     for j = 1:180,
         fprintf(fid,'%6.1f %6.1f %12.10e %12.10e %12.10e %12.10e %12.10e %12.10e\n', ...
         X1(j,i),Y1(j,i),Hrr(j,i),Hri(j,i),Htr(j,i),Hti(j,i),Hpr(j,i),Hpi(j,i));
     end;
 end;
         
 fclose all;

load c:\manoj\ocean\ECCO_Cord;
[X1,Y1] = meshgrid([-0.5;V_long],[90.5:-1:80.5 fliplr(U_lat') -80.5:-1:-90.5]);
[X2,Y2] = meshgrid(0:0.6:359.4,89.7:-.6:-89.7);
Hrr = [Hrr(1,:); Hrr; Hrr(end,:)];
Htr = [Htr(1,:); Htr; Htr(end,:)];
Hpr = [Hpr(1,:); Hpr; Hpr(end,:)];

Hrr_n = griddata(X1,Y1,[Hrr(:,end) Hrr],X2,Y2);
Htr_n = griddata(X1,Y1,[Htr(:,end) Htr],X2,Y2);
Hpr_n = griddata(X1,Y1,[Hpr(:,end) Hpr],X2,Y2);

[X2,Y2] = meshgrid(0:0.6:359.4,0.3:0.6:179.7);%stefan needed the data to be written with co-latitude

fid = fopen('c:\manoj\ocean\mf6\ECCO_2000_2006_MEAN_B_new.txt','wt');

fprintf(fid,'Colat   Long    Hr Ht Hp\n');
for j = 1:300,
for i = 1:600,
            fprintf(fid,'%10.1f %10.1f    %12.10e %12.10e %12.10e\n', ...
        Y2(j,i),X2(j,i),Hrr_n(j,i),Htr_n(j,i),Hpr_n(j,i));
    end;
end;
fclose all;
